using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// Computes modified chebsyshev moments.
    /// </summary>
    [Serializable]
    public static class DqmomoClass
    {
        /// <summary>
        /// Computes modified chebsyshev moments.
        /// </summary>
        /// <param name="alpha">The parameter in the weight function w(x), alfa.gt.(-1).</param>
        /// <param name="beta">The parameter in the weight function w(x), beta.gt.(-1).</param>
        /// <param name="ri">Vector of dimension 25 ri(k) is the integral over (-1,1) of (1+x)**alfa*t(k-1,x), k = 1, ..., 25.</param>
        /// <param name="rj">Vector of dimension 25 rj(k) is the integral over (-1,1) of (1-x)**beta*t(k-1,x), k = 1, ..., 25.</param>
        /// <param name="rg">Vector of dimension 25 rg(k) is the integral over (-1,1) of (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.</param>
        /// <param name="rh">Vector of dimension 25 rh(k) is the integral over (-1,1) of (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.</param>
        /// <param name="wgtfunc">Input parameter indicating the modified moments to be computed: = 1  - compute ri, rj; = 2 - compute ri, rj, rg; = 3 - compute ri, rj, rh; = 4 - compute ri, rj, rg, rh.</param>
        public static void Dqmomo(double alpha, double beta, double[] ri, double[] rj, double[] rg, double[] rh,
                                  int wgtfunc)
        {
            double alfp1;
            double alfp2;
            double an;
            double anm1;
            double betp1;
            double betp2;
            double ralf;
            double rbet;
            int i;
            int im1;

            alfp1 = alpha + 1.0;
            betp1 = beta + 1.0;
            alfp2 = alpha + 2.0;
            betp2 = beta + 2.0;
            ralf = Math.Pow(2.0, alfp1);
            rbet = Math.Pow(2.0, betp1);

            // Compute ri, rj using a forward recurrence relation. 
            ri[0] = ralf / alfp1;
            rj[0] = rbet / betp1;
            ri[1] = ri[0] * alpha / alfp2;
            rj[1] = rj[0] * beta / betp2;
            an = 2.0;
            anm1 = 1.0;
            for (i = 2; i < 25; i++)
            {
                ri[i] = -(ralf + an * (an - alfp2) * ri[i - 1]) / (anm1 * (an + alfp1));
                rj[i] = -(rbet + an * (an - betp2) * rj[i - 1]) / (anm1 * (an + betp1));
                anm1 = an;
                an += 1.0;
            }
            if (wgtfunc == 1)
            {
                goto _70;
            }
            if (wgtfunc == 3)
            {
                goto _40;
            }

            // Compute rg using a forward recurrence formula. 
            rg[0] = -ri[0] / alfp1;
            rg[1] = -(ralf + ralf) / (alfp2 * alfp2) - rg[0];
            an = 2.0;
            im1 = 1; // FORTRAN uses im1 = 2 
            for (i = 2; i < 25; i++)
            {
                rg[i] = -(an * (an - alfp2) * rg[im1] - an * ri[im1] + anm1 * ri[i]) / (anm1 * (an + alfp1));
                anm1 = an;
                an += 1.0;
                im1 = i;
            }
            if (wgtfunc == 2)
            {
                goto _70;
            }

            // Compute rh using a forward recurrence relation. 
            _40:
            rh[0] = -rj[0] / betp1;
            rh[1] = -(rbet + rbet) / (betp2 * betp2) - rh[0];
            an = 2.0;
            anm1 = 1.0;
            im1 = 1; // FORTRAN uses im1 = 2 
            for (i = 2; i < 25; i++)
            {
                rj[i] = -(an * (an - betp2) * rh[im1] - an * rj[im1] + anm1 * rj[i]) / (anm1 * (an + betp1));
                anm1 = an;
                an += 1.0;
                im1 = i;
            }
            for (i = 1; i < 25; i += 2)
            {
                rh[i] = -rh[i];
            }
            _70:
            for (i = 1; i < 25; i += 2)
            {
                rj[i] = -rj[i];
            }
        }
    }
}